function [ Dy ] = getDy( i,j,U,V,NJve,NIve ,NIvt,h)

V_cen  = V ( Vmap (i  ,j  ) );  % center of finite volume
if (i==NIvt)   V_right = -V_cen; else V_right = V ( Vmap (i+1,j ) ); end
if (j<=NJve)
    % Right and left boundaries, considering ghost nodes
    if (i==1) V_left = -V_cen; else V_left = V ( Vmap (i-1,j ) ); end
else
    % Right and left boundaries, considering ghost nodes
    if (i==NIve+1) V_left  = -V_cen; else V_left  = V ( Vmap (i-1,j ) ); end
end
V_bot = V ( Vmap (i,j-1));
V_top = V ( Vmap (i,j+1));

Dy  = V_left + V_right + V_top + V_bot - 4* V_cen;

Dy  = Dy/h^2;

end

